1 Πριν ξεκινήσουμε

Βεβαιώστε ότι το όνομα του project σας, όπως και το file path ΔΕΝ περιέχουν:
1. Ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Επί παραδείγματι το file path σας θα μπορούσε να είναι: C:/Environmental_Data_Analysis/Geographical_Data.

2 Φόρτωση βιβλιοθηκών

Ας φορτώσουμε όλες τις βιβλιοθήκες τις οποίες θα χρησιμοποιήσουμε σήμερα¹.

## ===========================================================================##
## Load the libraries
## ===========================================================================##
library(tidyverse)
library(raster)
library(sf)
library(ape)
library(phyloregion)
library(picante)
library(rtrees)
library(ggtree)
library(exactextractr)
library(caret)
library(psych)
library(PhyloMeasures)
library(phyloregion)
library(Matrix)
library(broom)
## ===========================================================================##

3 Σχέση έκτασης - αριθμού ειδών

Η σχέση αυτή θεωρείται ένας από τους ελάχιστους νόμους της οικολογίας, καθώς έχει ισχύ σε όλες τις ταξινομικές βαθμίδες και όλους τους οργανισμούς (Triantis et al., 2012). Με πολύ απλά λόγια, θεωρεί ότι όσο αυξάνεται η έκταση μιας περιοχής, τόσο περισσότερα taxa αυτή φιλοξενεί.

Όμως υπάρχουν και άλλοι παράγοντες (κλιματολογικοί, τοπογραφικοί, γεωλογικοί και εδαφικοί), οι οποίοι μπορούν να επηρεάσουν τα πρότυπα ποικιλότητας και κατανομής των ειδών στον χώρο και τον χρόνο.

4 Περιοχή μελέτης

Σήμερα θα ασχοληθούμε με τα Ryukyus Islands, τα οποία είναι περισσότερο γνωστά ως Nansei Islands (στα Ιαπωνικά: Nansei-shoto, δηλαδή “Νοτιοδυτικά Νησιά”), ανήκουν στην Ιαπωνία. Αποτελούν ένα νησιωτικό τόξο, το οποίο βρίσκεται σε γεωγραφικό πλάτος 24^ο - 31.0^οΝ και εκτείνεται σε μήκος 1300 km σε έναν άξονα ΝΔ-ΒΑ, ξεκινώντας από το Kyushu και καταλήγει στην Taiwan. Το νοτιότερο νησί είναι το Yonaguni, ενώ το μεγαλύτερο η Okinawa. Τα μεγαλύτερα σε μέγεθος νησιά είναι κυρίως ορεινά, ενώ τα μικρότερα είναι κοραλλιογενή. Καθώς εντοπίζεται στην ζώνη μετάβασης από το τροπικό στο θερμό εύκρατο κλίμα (δηλαδή, βρίσκεται στην υποτροπική ζώνη), οι κλιματικές συνθήκες είναι ήπιες στο αρχιπέλαγος αυτό, καθ’όλη την διάρκεια του έτους, με μέση θερμοκρασία 15^οC το χειμώνα και 28^οC το καλοκαίρι. Βέβαια, υπάρχουν διαβαθμίσεις: στο Β τμήμα του αρχιπελάγους, το κλίμα χαρακτηρίζεται ως υγρό υποτροπικό, ενώ στο Ν τμήμα του ως τροπικό (ταξινόμηση κατά Koppen). Εξαιτίας της γεωγραφικής τους θέσης, στα νησιά αυτά δεν υπάρχει ξηρή περίοδος και τα περισσότερα εξ αυτών δέχονται περισσότερα από 2000 m βροχής, λόγω των μουσώνων και των καλοκαιρινών τυφώνων. Ως εκ τούτου, στα Ryukyus Islands δεν εμφανίζεται η υποτροπική ξηρή ζώνη, η οποία διαχωρίζει τα εύκρατα δάση από τις υγρές τροπικές περιοχές σε άλλα μέρη του πλανήτη και η φυσική τους βλάστηση αποτελείται από πλατύφυλλα αειθαλή δαση (Good, 1845; Maekawa, 1974; Kira, 1991).

Γεωλογικά, το αρχιπέλαγος διαιρείται σε 3 τμήματα:
1. Northern Ryukyus (τα νησιά βορείως του Tokara Strait)
2. Central Ryukyus (τα νησιά μεταξύ του Tokara Strait και του Kerama Gap)
3. Southern Ryukyus (τα νησιά νοτίως του Kerama Gap)

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Τα Ryukyus είναι ηπειρωτικά νησιά και αποτελούν τις κορυφές της κορδιλλιέρας του αρχιπελάγους Ryukyus και έχουν δημιουργηθεί ως αποτέλεσμα της καταβύθισης της τεκτονικής πλάκας των Φιλιππίνων (Ota, 1998).

5 Εισαγωγή δεδομένων

Θα χρησιμοποιήσουμε τα ελεύθερα προσβάσιμα δεδομένα από τους Nakamura et al. (2009) για ορισμένα από τα νησιά αυτά. Τώρα ας εισαγάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε:

## ===========================================================================##
## Load our data
## ===========================================================================##
plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", 
    row.names = 1)
## ===========================================================================##

Ας δούμε πόσες στήλες και σειρές έχει:

## ===========================================================================##
## Check their dimensions
## ===========================================================================##
dim(plants)

## [1]   26 1815

Και τι είδους δεδομένα περιέχει:

## ===========================================================================##
## Check their dimensions
## ===========================================================================##
str(plants)

## 'data.frame':    26 obs. of  1815 variables:
##  $ Codonacanthus_pauciflorus               : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Dicliptera_chinensis                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hemigraphis_reptans                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hygrophila_salicifolia                  : int  0 0 0 1 0 0 0 1 0 0 ...
##  $ Justicia_procumbens                     : int  0 0 1 1 1 1 1 1 0 1 ...
##  $ Lepidagathis_formosensis                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Lepidagathis_inaequalis                 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Staurogyne_concinnula                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Strobilanthes_glandulifer               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Strobilanthes_tashiroi                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_capillipes                         : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Acer_insulare                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_oblongum                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_rufinerve                          : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Acer_sieboldianum                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_hypoleuca                     : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_polygama                      : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_rufa                          : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Saurauria_roxburghii                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Mollugo_pentaphylla                     : int  0 0 1 1 1 0 0 0 0 1 ...
##  $ Tetragonia_tetragonioides               : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ Alangium_premnifolium                   : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Alisma_canaliculatum                    : int  0 0 0 1 0 0 0 1 0 1 ...
##  $ Caldesia_reniformis                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Sagittaria_aginashi                     : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Sagittaria_pygmaea                      : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Sagittaria_trifolia                     : int  0 0 0 1 1 1 0 0 0 0 ...
##  $ Achyranthes_aspera                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Achyranthes_bidentata                   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Achyranthes_longifolia                  : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Deeringia_polysperma                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Philoxerus_wrightii                     : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Crinum_asiaticum                        : int  1 0 1 1 1 1 1 1 1 1 ...
##  $ Lycoris_traubii                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_crenata                            : int  0 0 1 1 1 0 1 1 0 0 ...
##  $ Rhus_ambigua                            : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Rhus_javanica                           : int  0 1 1 1 1 1 0 0 0 0 ...
##  $ Rhus_succedanea                         : int  0 0 1 1 1 1 0 1 0 0 ...
##  $ Polyalthia_liukiuensis                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Anodendron_affine                       : int  0 1 1 1 1 1 0 0 0 0 ...
##  $ Cerbera_manghas                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ecdysanthera_utilis                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Neisosperma_oppositifolia               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Parsonsia_laevigata                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Trachelospermum_asiaticum               : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ Trachelospermum_jasminoides             : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Ilex_dimorphophylla                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_ficoidea                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_goshiensis                         : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_hayataiana                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_integra                            : int  1 1 1 1 1 1 1 1 1 0 ...
##  $ Ilex_kusanoi                            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_liukiuensis                        : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_macrocarpa                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_macropoda                          : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Ilex_maximowicziana                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_pedunculosa                        : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_purpurea                           : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Ilex_rotunda                            : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Acorus_gramineus                        : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Alocasia_cucullata                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alocasia_gageana                        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alocasia_odora                          : int  0 0 1 1 1 1 1 1 0 1 ...
##  $ Amorphophallus_kiusianus                : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Arisaema_heterocephalum                 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Arisaema_japonicum                      : int  0 0 1 1 1 0 0 1 0 0 ...
##  $ Arisaema_kawashimae                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Arisaema_longipedunculatum              : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_ringens                        : int  0 0 1 1 1 1 1 1 1 1 ...
##  $ Arisaema_robustum                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_sazensoo                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_thunbergii                     : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Arisaema_yakualpinum                    : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Pinellia_ternata                        : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Pinellia_tripartia                      : int  0 0 0 0 1 0 0 0 1 1 ...
##  $ Rhaphidophora_liukiuensis               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Rhaphidophora_pinnata                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Typhonium_divaricatum                   : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Aralia_elata                            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Dendropanax_trifidus                    : int  0 1 1 1 1 1 1 1 1 0 ...
##  $ Fatsia_japonica                         : int  1 1 1 1 1 1 1 1 1 0 ...
##  $ Hedera_rhombea                          : int  0 0 1 1 1 0 1 1 0 0 ...
##  $ Kalopanax_pictus                        : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Kalopanax_septemlobus                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Schefflera_octophylla                   : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Aristolochia_debilis                    : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Aristolochia_kaempferi                  : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Aristolochia_liukiuensis                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Aristolochia_zollingeriana              : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_caudigerum                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_celsum                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_dissitum                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_fudsinoi                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_gelasinum                        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_gusuk                            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_hatsushimae                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_hirsuitsepalum                   : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Asarum_lutchuense                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_monodoriflorum                   : int  0 0 0 0 0 0 0 0 0 0 ...
##   [list output truncated]

6 Μετατροπή μήτρας παρουσίας/απουσίας

Όπως θα διαπιστώσατε, τα δεδομένα αυτά αποτελούν μια μήτρα παρουσίας/απουσίας ειδών (presence/absence matrix). Τέτοιου είδους δεδομένα είναι πολύ χρήσιμα για αναλύσεις σχετικά με την β-ποικιλότητα, όμως σήμερα εμείς χρειαζόμαστε μόνο τον αριθμό των ειδών σε ένα συγκεκριμένο νησί και όχι ποια είδη βρίσκονται σε κάθε νησί.

Ένας τρόπος προκειμένου να μετατρέψουμε μια μήτρα παρουσίας/απουσία ειδών στον συνολικό αριθμό ειδών τα οποία απαντώνται σε κάθε νησί, είναι ο ακόλουθος:

## ===========================================================================##
## Convert the P/A matrix to SR
## ===========================================================================##
ryukyus <- rowSums(plants) %>% enframe() %>% dplyr::rename(Island = name, SR = value)

ryukyus

7 Εισαγωγή γεωχωρικών δεδομένων

Καθώς έχουμε πλέον στη μορφή που θέλουμε τα βιοτικά μας δεδομένα, επόμενο μας βήμα είναι να συγκεντρώσουμε αβιοτικά δεδομένα για τα νησιά αυτά.

## ===========================================================================##
## Create a vector containing just the file paths
## ===========================================================================##
shps <- fs::dir_ls("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/", 
    glob = "*.shp")
## ===========================================================================##


## ===========================================================================##
## Create a vector containing just the island names
## ===========================================================================##
island_nms <- shps %>% str_replace_all("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/", 
    "") %>% str_replace_all("\\.shp", "")
## ===========================================================================##


## ===========================================================================##
## Bulk load the shapefiles
## ===========================================================================##
ryukyus_shps <- shps %>% purrr::map(read_sf) %>% lapply(., dplyr::select, NAME_1, 
    geometry) %>% do.call(rbind, .) %>% mutate(Island = island_nms) %>% dplyr::select(Island, 
    geometry)
## ===========================================================================##

## ===========================================================================##
## Load, crop and mask the rasters
## ===========================================================================##
preds <- stack(stack(list.files(path = "I:/SDM DATA/WORLDCLIM/Current/bio_2-5m_bil/", 
    pattern = ".bil", recursive = T, full.names = T)) %>% crop(., ryukyus_shps, snap = "in") %>% 
    mask(., ryukyus_shps), stack(list.files(path = "I:/SDM DATA/Envirem_Predictors/Eurasia_current_2.5arcmin_geotiff/", 
    pattern = ".tif", recursive = T, full.names = T)) %>% crop(., ryukyus_shps, snap = "in") %>% 
    mask(., ryukyus_shps))
## ===========================================================================##

## ===========================================================================##
## Load, crop and mask the altitudinal raster
## ===========================================================================##
altitude <- raster("I:/SDM DATA/Altitude/elevation_1KMmd_GMTEDmd.tif") %>% crop(., 
    ryukyus_shps, snap = "in") %>% mask(., ryukyus_shps)
## ===========================================================================##

8 Εξαγωγή γεωχωρικών δεδομένων

Τώρα μπορούμε να εξάγουμε τα αβιοτικά αυτά δεδομένα για κάθε νησί που εμπεριέχεται στην ανάλυση μας:

## ===========================================================================##
## Extract abiotic info
## ===========================================================================##
preds_ryukyus <- exact_extract(preds, ryukyus_shps, "median") %>% mutate(Island = ryukyus_shps$Island) %>% 
    full_join(., exact_extract(altitude, ryukyus_shps, "median") %>% enframe() %>% 
        mutate(Island = ryukyus_shps$Island) %>% rename(Elevation = value) %>% dplyr::select(Island, 
        Elevation)) %>% dplyr::select(Island, everything()) %>% as_tibble()

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%

## ============================================================================##
## Give them proper names
## ============================================================================##
col_names <- preds_ryukyus %>% colnames() %>% enframe() %>% mutate(value = str_replace_all(value, 
    "median\\.", "")) %>% mutate(value = str_replace_all(value, "current_2\\.5arcmin_", 
    "")) %>% pull(value)
## ===========================================================================##


## ============================================================================##
## Overwrite names and remove Row 1
## ============================================================================##
vars <- preds_ryukyus %>% set_names(col_names) %>% dplyr::select(Island, everything()) %>% 
    as_tibble()

9 Αντιμετώπιση προβλημάτων

Εάν παρατηρήσετε τα δεδομένα για τα νησιά που μας ενδιαφέρουν, θα δείτε², ότι για ορισμένα νησιά δεν έχουμε καθόλου (NA) δεδομένα. Αυτό μπορεί να οφείλεται είτε στο γεγονός ότι τα νησιά αυτά μπορεί να είναι πολύ μικρά, οπότε να μην υπάρχουν διαθέσιμα δεδομένα στο μέγεθος των ψηφιδοπλεγμάτων τα οποία χρησιμοποιούμε (resolution - ακρίβεια), είτε εξαιτίας του γεγονότος ότι χρησιμοποιούμε σχετικά αδρά ψηφιδοπλέγματα να μην έχουμε δεδομένα για αυτά τα νησιά, ανεξαρτήτως του μεγέθους τους³.

Υπάρχουν δύο αποδεκτές λύσεις:

• Μετασχηματισμός των δεδομένων μας
  
• Αφαίρεση των προβληματικών δεδομένων

Και οι δύο προσεγγίσεις έχουν τα πλεονεκτήματα και τα μειονεκτήματα τους. Η πρώτη προσέγγιση είναι προτιμότερη όταν έχουμε πλήθος διαθέσιμων δεδομένων (>5000 δείγματα) και το ποσοστό των προβληματικών δεδομένων είναι μικρό (συνήθως < 1-2%). Η δεύτερη προσέγγιση θεωρείται γενικά ασφαλέστερη, διότι δεν δημιουργούμε εκ του μηδενός δεδομένα. Σήμερα θα ακολουθήσουμε την δεύτερη προσέγγιση, όμως θα δούμε πως μπορούμε να μετασχηματίσουμε τα δεδομένα μας.

## ============================================================================##
## Transform the missing data
## ============================================================================##
imputed_vars <- predict(preProcess(as.data.frame(vars[, 2:36]), method = c("medianImpute")), 
    vars)

## ============================================================================##
## First drop any rows with NA values
## ============================================================================##
drop_vars <- vars %>% drop_na()
## ============================================================================##


## ============================================================================##
## Then check which islands have NA values in any of the variables
## ============================================================================##
setdiff(imputed_vars$Island, drop_vars$Island)

##  [1] "Akusekijima"  "Gajajima"     "Ioujima"      "Kodakarajima" "Kuchinoshima"
##  [6] "Tairajima"    "Takarajima"   "Takeshima"    "Yokoatejima"  "Yonagunijima"

## ============================================================================##
## Create a symbol to filter out all the Islands having NA values
## ============================================================================##
"%!in%" <- function(x, y) !(x %in% y)
## ============================================================================##


## ============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
## ============================================================================##
plants_new <- plants %>% mutate(Island = rownames(.)) %>% as_tibble() %>% dplyr::select(Island, 
    everything()) %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>% 
    filter(Island %!in% setdiff(imputed_vars$Island, drop_vars$Island)) %>% pivot_wider(names_from = Taxon, 
    values_from = PA)
## ============================================================================##


## ============================================================================##
## See which taxa are left after removing the problematic islands
## ============================================================================##
keepers <- colSums(plants_new[, -1]) %>% enframe() %>% filter(value > 0)
## ============================================================================##


## ============================================================================##
## Finally create the new PA matrix
## ============================================================================##
plants_new_pa <- plants_new %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>% 
    filter(Taxon %in% keepers$name) %>% pivot_wider(names_from = Taxon, values_from = PA)
## ============================================================================##

10 Κατασκευή φυλογενετικού δένδρου

Προκειμένου να μπορέσουμε να υπολογίσουμε την φυλογενετική α-ποικιλότητα για το αρχιπέλαγος της Ιαπωνίας, θα χρειαστεί να κατασκευάσουμε πρώτα το φυλογενετικό δένδρο των ειδών που απαρτίζουν τις νησιωτικές φυτοκοινότητες. Για τον σκοπό αυτό, θα χρησιμοποιήσουμε εντολές από τις βιβλιοθήκες tidyr, dplyr & rtrees. Η τελευταία βιβλιοθήκη περιλαμβάνει τα φυλογενετικά δένδρα για αρκετούς φυτικούς οργανισμούς, καθώς και για τους ιχθύες, τα πτηνά και τα θηλαστικά. Θα χρησιμοποιήσουμε μόνο τα taxa τα οποία απαντώνται στα 16 νησιά του αρχιπελάγους της Ιαπωνίας που δεν είχαν κάποιο πρόβλημα με τις αβιοτικές μεταβλητές. Τα taxa αυτά βρίσκονται στον πίνακα που ονομάσαμε keepers. Θα χρειαστεί να προσθέσουμε στον πίνακα αυτόν το όνομα του γένους για κάθε φυτικό taxon το οποίο απαντάται στα 16 αυτά νησιά των Ryukyus και στη συνέχεια, θα χρησιμοποιήσουμε αυτόν τον πίνακα για να κατασκευάσουμε το φυλογενετικό μας δένδρο. Ακόμα και εάν το φυλογενετικό μας δένδρο περιέχει πολυτομίες, είναι ασφαλές να το χρησιμοποιήσουμε για μακροοικολογικές αναλύσεις (Li et al., 2019).

Αναλόγως των δυνατοτήτων του Η/Υ σας, υπάρχει περίπτωση η ακόλουθη εντολή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή: jpn_tree <- readRDS('Data/RDS/16 Ryukyus tree.rds')

##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- keepers %>% 
  separate(name, into = c('genus', 'binomial'), sep = '_', remove = F) %>% 
  rename(species = name) %>% 
  dplyr::select(species, genus) %>% 
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)
##============================================================================##

Στη συνέχεια, μπορούμε να αναπαραστήσουμε γραφικά το φυλογενετικό δένδρο, δίνοντας διαφορετικό χρώμα σε κάθε γένος:

## ============================================================================##
## Plot the tree
## ============================================================================##
groupInfo <- split(jpn_tree$tip.label, gsub("_\\w+", "", jpn_tree$tip.label))

groups_jpn <- groupOTU(jpn_tree, groupInfo)
ggtree(groups_jpn, aes(color = group), layout = "circular") + geom_tiplab2(size = 0.5, 
    align = T)

## ============================================================================##
## Combine the trees with the P/A data
## ============================================================================##

# check for mismatches/missing species
combined_japan <- picante::match.phylo.comm(jpn_tree, plants_new_pa[, 2:1800])

## [1] "Dropping taxa from the community because they are not present in the phylogeny:"
##  [1] "Saurauria_roxburghii"       "Carmona_retusa"            
##  [3] "Saionia_shinzatoi"          "Crataeva_falcata"          
##  [5] "Pseudostellaris_heterantha" "Dispanthus_uniflorus"      
##  [7] "Siegesbeckia_glabrescens"   "Siegesbeckia_orientalis"   
##  [9] "Cardamina_flexuosa"         "Pogonantherum_crinitum"    
## [11] "Spodipogon_cotulifer"       "Eusteralis_stellata"       
## [13] "Perillura_reptans"          "Akebi_quinata"             
## [15] "Dumbaria_villosa"           "Alectrorurus_yedoensis"    
## [17] "Caphalanthera_falcata"      "Heterozeuxine_nervosa"     
## [19] "Heterozeuxine_odorata"      "Luisaerides_liukiuensis"   
## [21] "Vexillabium_yakushimense"   "Hydyotis_verticillata"     
## [23] "phtheirospermum_japonicum"  "Cnidum_japonicum"          
## [1] "Dropping tips from the tree because they are not present in the community data:"
## [1] "Phtheirospermum_japonicum"

## ============================================================================##
## The resulting object is a list with $phy and $comm elements replace our
## original data with the sorted/matched data
## ============================================================================##
phy_japan <- combined_japan$phy

comm_japan <- combined_japan$comm
## ============================================================================##

## ============================================================================##
## Compute phylogenetic alpha diversity
## ============================================================================##

## For the estimation of SES values
japan_pd <- pd.query(phy_japan, comm_japan, TRUE)

## For the p-values
japan_pd_p <- pd.pvalues(phy_japan, comm_japan, reps = 1000)
## ============================================================================##

11 Τελικό dataset

Τώρα πρέπει να κρατήσουμε μόνο εκείνα τα νησιά που δεν έχουν προβληματικά δεδομένα, να τα ενώσουμε με τα αβιοτικά και τα φυλογενετικά δεδομένα, καθώς και να υπολογίσουμε την έκταση του εκάστοτε νησιού:

##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>% 
  filter(Island %in% plants_new_pa$Island) %>% 
  full_join(., drop_vars) %>% 
  full_join(., ryukyus_shps) %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  st_as_sf() %>% 
  
  ## Lambert Azimuthal Equal Area projection
  st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>% 
  
  ## Create the Area and PD variables
  mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
         PD = japan_pd,
         PD_P = japan_pd_p) %>% 
  
  ## Transform back to WGS84
  st_transform(4326) %>% 
  
  dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>% 
  
  ## It has zero-variance
  dplyr::select(-aridityIndexThornthwaite) %>% 
  
  ## Drop the geometry
  st_drop_geometry() %>% 
  
  ## Join these data with the Maximum and Minimum distance between each island
  ## pair
  inner_join(., 
             
             ## Calculate the Maximum distance between each island pair
             apply((st_distance(ryukyus_shps %>% 
                                     st_centroid(), 
                                   ryukyus_shps %>% 
                                     st_centroid())/1e+3) %>% 
                        as_tibble() %>% 
                        set_names(ryukyus_shps$Island) %>% 
                        units::drop_units() %>% 
                        na_if(., 0), 1, max, na.rm = T) %>% 
               enframe() %>% 
               rename(MaxD = value) %>% 
               mutate(Island = ryukyus_shps$Island) %>% 
               dplyr::select(Island, MaxD) %>% 
               
               full_join(., 
                         
                         ## Calculate the Minimum distance between each island
                         ## pair
                         apply((st_distance(ryukyus_shps %>% 
                                                 st_centroid(), 
                                               ryukyus_shps %>% 
                                                 st_centroid())/1e+3) %>% 
                                    as_tibble() %>% 
                                    set_names(ryukyus_shps$Island) %>% 
                                    units::drop_units() %>% 
                                    na_if(., 0), 1, min, na.rm = T) %>% 
                           enframe() %>% 
                           rename(MinD = value) %>% 
                           mutate(Island = ryukyus_shps$Island) %>% 
                           dplyr::select(Island, MinD)))
##============================================================================##

ryukyus_new

12 Προκαταρτικές αναλύσεις

## ============================================================================##
## Normality check
## ============================================================================##
normality <- lapply(ryukyus_new[, 2:42], shapiro.test)
## ============================================================================##

normality

## $SR
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93472, p-value = 0.2893
## 
## 
## $PD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.96192, p-value = 0.6968
## 
## 
## $PD_P
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.77558, p-value = 0.001314
## 
## 
## $Area
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.75456, p-value = 0.0007212
## 
## 
## $Elevation
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.72274, p-value = 0.0003039
## 
## 
## $bio1
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93352, p-value = 0.2771
## 
## 
## $bio10
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.85281, p-value = 0.01499
## 
## 
## $bio11
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93575, p-value = 0.3001
## 
## 
## $bio12
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.90384, p-value = 0.09261
## 
## 
## $bio13
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.86312, p-value = 0.02138
## 
## 
## $bio14
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92022, p-value = 0.1701
## 
## 
## $bio15
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92241, p-value = 0.1845
## 
## 
## $bio16
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91524, p-value = 0.1414
## 
## 
## $bio17
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88965, p-value = 0.055
## 
## 
## $bio18
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92892, p-value = 0.2344
## 
## 
## $bio19
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.89273, p-value = 0.06153
## 
## 
## $bio2
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94329, p-value = 0.3913
## 
## 
## $bio3
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91761, p-value = 0.1544
## 
## 
## $bio4
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94292, p-value = 0.3863
## 
## 
## $bio5
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.82406, p-value = 0.005789
## 
## 
## $bio6
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93515, p-value = 0.2938
## 
## 
## $bio7
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94416, p-value = 0.4031
## 
## 
## $bio8
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88975, p-value = 0.0552
## 
## 
## $bio9
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92902, p-value = 0.2352
## 
## 
## $annualPET
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88922, p-value = 0.05414
## 
## 
## $climaticMoistureIndex
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92624, p-value = 0.2124
## 
## 
## $continentality
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93748, p-value = 0.3193
## 
## 
## $embergerQ
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91509, p-value = 0.1406
## 
## 
## $growingDegDays0
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92894, p-value = 0.2345
## 
## 
## $growingDegDays5
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92894, p-value = 0.2345
## 
## 
## $maxTempColdest
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92989, p-value = 0.2428
## 
## 
## $minTempWarmest
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88613, p-value = 0.04841
## 
## 
## $monthCountByTemp10
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.46956, p-value = 1.199e-06
## 
## 
## $PETColdestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9442, p-value = 0.4037
## 
## 
## $PETDriestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92498, p-value = 0.2028
## 
## 
## $PETseasonality
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92548, p-value = 0.2065
## 
## 
## $PETWarmestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.87995, p-value = 0.03873
## 
## 
## $PETWettestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91681, p-value = 0.1498
## 
## 
## $thermicityIndex
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93209, p-value = 0.2631
## 
## 
## $MaxD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91265, p-value = 0.1284
## 
## 
## $MinD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8222, p-value = 0.005454

## ============================================================================##
## Correlation check
## ============================================================================##
corr.test(ryukyus_new[, 5:42], method = "spearman")

## Call:corr.test(x = ryukyus_new[, 5:42], method = "spearman")
## Correlation matrix 
##                        Area Elevation  bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area                   1.00      0.06  0.17  0.11  0.14 -0.01 -0.22  0.01 -0.32
##                       bio16 bio17 bio18 bio19  bio2  bio3  bio4  bio5  bio6
## Area                  -0.15 -0.07  0.25 -0.05 -0.22  0.11 -0.22  0.10  0.10
##                        bio7  bio8  bio9 annualPET climaticMoistureIndex
## Area                  -0.20  0.24  0.00      0.11                 -0.01
##                       continentality embergerQ growingDegDays0 growingDegDays5
## Area                           -0.22      0.32            0.16            0.16
##                       maxTempColdest minTempWarmest monthCountByTemp10
## Area                            0.15           0.09              -0.09
##                       PETColdestQuarter PETDriestQuarter PETseasonality
## Area                               0.18            -0.05          -0.44
##                       PETWarmestQuarter PETWettestQuarter thermicityIndex  MaxD
## Area                              -0.06             -0.23            0.14  0.04
##                        MinD
## Area                   0.49
##  [ reached getOption("max.print") -- omitted 37 rows ]
## Sample Size 
## [1] 16
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##                       Area Elevation bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area                  0.00      1.00 1.00  1.00  1.00  1.00  1.00  1.00  1.00
##                       bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7
## Area                   1.00  1.00  1.00  1.00 1.00 1.00 1.00 1.00 1.00 1.00
##                       bio8 bio9 annualPET climaticMoistureIndex continentality
## Area                  1.00 1.00      1.00                  1.00           1.00
##                       embergerQ growingDegDays0 growingDegDays5 maxTempColdest
## Area                       1.00            1.00            1.00           1.00
##                       minTempWarmest monthCountByTemp10 PETColdestQuarter
## Area                            1.00               1.00              1.00
##                       PETDriestQuarter PETseasonality PETWarmestQuarter
## Area                              1.00           1.00              1.00
##                       PETWettestQuarter thermicityIndex MaxD MinD
## Area                               1.00            1.00 1.00 1.00
##  [ reached getOption("max.print") -- omitted 37 rows ]
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

## ============================================================================##
## Correlation check
## ============================================================================##
unconc <- usdm::vifcor(ryukyus_new %>% as.data.frame() %>% dplyr::select(Area:MinD), 
    th = 0.7)
unconc <- unconc@results$Variables[unconc@results$VIF <= 5]
unconc <- droplevels(factor(unconc))
## ============================================================================##

## ============================================================================##
## Visualization
## ============================================================================##
car::scatterplotMatrix(ryukyus_new[, 5:10], spread = F, smoother.args = list(lty = 2), 
    main = "Scatter Plot Matrix")

## ============================================================================##
## Log-transform the data
## ============================================================================##
ryukyus_log <- ryukyus_new %>% dplyr::select(Area, bio14, bio8, monthCountByTemp10, 
    MaxD, MinD) %>% mutate(SR = ryukyus_new$SR) %>% log10(.) %>% mutate(Island = ryukyus_new$Island, 
    PD = ryukyus_new$PD) %>% dplyr::select(Island, SR, PD, everything())
## ============================================================================##

13 Κυρίως ανάλυση

## ============================================================================##
## ISAR models
## ============================================================================##
model1 <- lm(SR ~ Area, data = ryukyus_log)
model2 <- lm(PD ~ Area, data = ryukyus_log)
## ============================================================================##

## ============================================================================##
## Check the summary
## ============================================================================##
summary(model1)

## 
## Call:
## lm(formula = SR ~ Area, data = ryukyus_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.169622 -0.049189 -0.001934  0.072756  0.107531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.17927    0.08117  26.848 1.93e-13 ***
## Area         0.28263    0.03750   7.537 2.72e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08492 on 14 degrees of freedom
## Multiple R-squared:  0.8023, Adjusted R-squared:  0.7882 
## F-statistic: 56.81 on 1 and 14 DF,  p-value: 2.722e-06

## ============================================================================##
## Check the summary
## ============================================================================##
summary(model2)

## 
## Call:
## lm(formula = PD ~ Area, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02464 -0.50739  0.04005  0.60604  1.37789 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.4217     0.6971  -3.474  0.00372 **
## Area          0.5013     0.3221   1.557  0.14188   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7293 on 14 degrees of freedom
## Multiple R-squared:  0.1475, Adjusted R-squared:  0.08665 
## F-statistic: 2.423 on 1 and 14 DF,  p-value: 0.1419

## ============================================================================##
## Create the residuals
## ============================================================================##
res_sr <- rstandard(model1)
res_pd <- rstandard(model2)
## ============================================================================##

## ============================================================================##
## Check the estimate and the p-value
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>% full_join(., lm(PD ~ 
    Area, data = ryukyus_log) %>% tidy()) %>% mutate(Metric = rep(c("SR", "PD"), 
    each = 2)) %>% dplyr::select(Metric, everything())
## ============================================================================##

## ============================================================================##
## Check the metrics
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>% full_join(., lm(PD ~ 
    Area, data = ryukyus_log) %>% glance()) %>% mutate(Metric = c("SR", "PD")) %>% 
    dplyr::select(Metric, everything())
## ============================================================================##


## ============================================================================##
## Check the residuals
## ============================================================================##
resids <- lm(SR ~ Area, data = ryukyus_log) %>% augment() %>% mutate(Island = ryukyus_log$Island) %>% 
    rename(SR_res = .std.resid) %>% full_join(., lm(PD ~ Area, data = ryukyus_log) %>% 
    augment() %>% mutate(Island = ryukyus_log$Island) %>% rename(PD_res = .std.resid), 
    by = "Island") %>% dplyr::select(Island, SR, SR_res, PD, PD_res)
## ============================================================================##

## ============================================================================##
## Plot
## ============================================================================##
ggplot(ryukyus_log, aes(res_sr, res_pd)) + geom_point() + geom_text(aes(label = Island), 
    fontface = "bold", size = 2.8, vjust = -1.2, hjust = 0.5) + geom_hline(yintercept = 0) + 
    geom_vline(xintercept = 0) + theme(legend.position = "none") + xlab("SR") + ylab("PD")

## ============================================================================##
## Quick save
## ============================================================================##
ggsave("SR-PD.png", width = 20, height = 20, units = "cm", dpi = 600)

13.2 Μοντελοποίηση - Πολλαπλή γραμμική παλινδρόμηση

13.2.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter

Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 2 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για την ταξινομική και ένα για την φυλογενετική ποικιλότητα σε σχέση με όλες τις ανεξάρτητες, μη συγγραμικές μεταβλητές μας. Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:

## ============================================================================##
## The enter method
## ============================================================================##
enter1 <- lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)

enter2 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)
## ============================================================================##

ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας.

Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:

## ============================================================================##
## Check the summary
## ============================================================================##
summary(enter1)

## 
## Call:
## lm(formula = SR ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MaxD + MinD, data = ryukyus_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.135024 -0.049051  0.003223  0.031492  0.114373 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.09862    2.01050   0.546   0.5980    
## Area                0.27193    0.04050   6.715  8.7e-05 ***
## bio14               0.17172    0.48562   0.354   0.7318    
## bio8                1.04033    0.59001   1.763   0.1117    
## monthCountByTemp10 -1.52961    0.81899  -1.868   0.0946 .  
## MaxD                0.02294    0.30523   0.075   0.9417    
## MinD               -0.09494    0.13400  -0.708   0.4966    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08294 on 9 degrees of freedom
## Multiple R-squared:  0.8788, Adjusted R-squared:  0.7979 
## F-statistic: 10.87 on 6 and 9 DF,  p-value: 0.001092

summary(enter2)

## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MaxD + MinD, data = ryukyus_log)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8957 -0.2648 -0.0292  0.4049  0.7872 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -20.5226    14.1621  -1.449   0.1812  
## Area                 0.7156     0.2853   2.509   0.0334 *
## bio14                7.5665     3.4207   2.212   0.0543 .
## bio8                 5.3023     4.1561   1.276   0.2340  
## monthCountByTemp10  -5.9614     5.7691  -1.033   0.3284  
## MaxD                -0.3244     2.1501  -0.151   0.8834  
## MinD                -2.1573     0.9439  -2.285   0.0481 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5842 on 9 degrees of freedom
## Multiple R-squared:  0.6484, Adjusted R-squared:  0.4139 
## F-statistic: 2.766 on 6 and 9 DF,  p-value: 0.08269

13.2.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης

Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών

## ============================================================================##
## The stepAIC method
## ============================================================================##
stepm1 <- MASS::stepAIC(lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + 
    MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL, 
    steps = 1000, k = 2)

stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + 
    MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL, 
    steps = 1000, k = 2)

summary(stepm1)

Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της φυλογενετικής ποικιλότητας;

ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared

## ============================================================================##
## The solution
## ============================================================================##
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + 
    MinD, data = ryukyus_log), scale = 0, direction = "both", trace = 1, keep = NULL, 
    steps = 1000, k = 2)

## Start:  AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - MaxD                1   0.00777 3.0795 -14.3653
## - monthCountByTemp10  1   0.36444 3.4361 -12.6119
## <none>                            3.0717 -12.4058
## - bio8                1   0.55550 3.6272 -11.7461
## - bio14               1   1.66988 4.7416  -7.4595
## - MinD                1   1.78274 4.8544  -7.0831
## - Area                1   2.14766 5.2194  -5.9234
## 
## Step:  AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - monthCountByTemp10  1   0.38000 3.4595 -14.5036
## <none>                            3.0795 -14.3653
## - bio8                1   0.56038 3.6398 -13.6903
## + MaxD                1   0.00777 3.0717 -12.4058
## - bio14               1   1.73902 4.8185  -9.2021
## - MinD                1   1.85944 4.9389  -8.8071
## - Area                1   2.13989 5.2194  -7.9234
## 
## Step:  AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - bio8                1   0.33274 3.7922 -15.0342
## <none>                            3.4595 -14.5036
## + monthCountByTemp10  1   0.38000 3.0795 -14.3653
## + MaxD                1   0.02333 3.4361 -12.6119
## - bio14               1   1.46668 4.9261 -10.8485
## - MinD                1   1.81014 5.2696  -9.7701
## - Area                1   2.38991 5.8494  -8.1001
## 
## Step:  AIC=-15.03
## PD ~ Area + bio14 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## <none>                            3.7922 -15.0342
## + bio8                1   0.33274 3.4595 -14.5036
## + monthCountByTemp10  1   0.15236 3.6398 -13.6903
## + MaxD                1   0.00319 3.7890 -13.0477
## - bio14               1   1.21835 5.0106 -12.5767
## - MinD                1   1.62428 5.4165 -11.3303
## - Area                1   2.38796 6.1802  -9.2199

summary(stepm2)

## 
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99647 -0.28430 -0.05602  0.32450  0.81760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -12.8387     6.5774  -1.952   0.0747 .
## Area          0.7480     0.2721   2.749   0.0176 *
## bio14         5.8728     2.9910   1.963   0.0732 .
## MinD         -1.5102     0.6661  -2.267   0.0427 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared:  0.5659, Adjusted R-squared:  0.4573 
## F-statistic: 5.214 on 3 and 12 DF,  p-value: 0.01554

## Start:  AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - MaxD                1   0.00777 3.0795 -14.3653
## - monthCountByTemp10  1   0.36444 3.4361 -12.6119
## <none>                            3.0717 -12.4058
## - bio8                1   0.55550 3.6272 -11.7461
## - bio14               1   1.66988 4.7416  -7.4595
## - MinD                1   1.78274 4.8544  -7.0831
## - Area                1   2.14766 5.2194  -5.9234
## 
## Step:  AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - monthCountByTemp10  1   0.38000 3.4595 -14.5036
## <none>                            3.0795 -14.3653
## - bio8                1   0.56038 3.6398 -13.6903
## + MaxD                1   0.00777 3.0717 -12.4058
## - bio14               1   1.73902 4.8185  -9.2021
## - MinD                1   1.85944 4.9389  -8.8071
## - Area                1   2.13989 5.2194  -7.9234
## 
## Step:  AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - bio8                1   0.33274 3.7922 -15.0342
## <none>                            3.4595 -14.5036
## + monthCountByTemp10  1   0.38000 3.0795 -14.3653
## + MaxD                1   0.02333 3.4361 -12.6119
## - bio14               1   1.46668 4.9261 -10.8485
## - MinD                1   1.81014 5.2696  -9.7701
## - Area                1   2.38991 5.8494  -8.1001
## 
## Step:  AIC=-15.03
## PD ~ Area + bio14 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## <none>                            3.7922 -15.0342
## + bio8                1   0.33274 3.4595 -14.5036
## + monthCountByTemp10  1   0.15236 3.6398 -13.6903
## + MaxD                1   0.00319 3.7890 -13.0477
## - bio14               1   1.21835 5.0106 -12.5767
## - MinD                1   1.62428 5.4165 -11.3303
## - Area                1   2.38796 6.1802  -9.2199

## 
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99647 -0.28430 -0.05602  0.32450  0.81760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -12.8387     6.5774  -1.952   0.0747 .
## Area          0.7480     0.2721   2.749   0.0176 *
## bio14         5.8728     2.9910   1.963   0.0732 .
## MinD         -1.5102     0.6661  -2.267   0.0427 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared:  0.5659, Adjusted R-squared:  0.4573 
## F-statistic: 5.214 on 3 and 12 DF,  p-value: 0.01554

13.2.2.1 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

## ============================================================================##
## Compare the two models
## ============================================================================##
AIC(enter1, stepm1)

Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;

13.2.2.2 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών

## ============================================================================##
## A function for the relative importance
## ============================================================================##
relweights <- function(fit, ...) {
    R <- cor(fit$model)
    nvar <- ncol(R)
    rxx <- R[2:nvar, 2:nvar]
    rxy <- R[2:nvar, 1]
    svd <- eigen(rxx)
    evec <- svd$vectors
    ev <- svd$values
    delta <- diag(sqrt(ev))
    lambda <- evec %*% delta %*% t(evec)
    lambdasq <- lambda^2
    beta <- solve(lambda) %*% rxy
    rsquare <- colSums(beta^2)
    rawwgt <- lambdasq %*% beta^2
    import <- (rawwgt/rsquare) * 100
    import <- as.data.frame(import)
    row.names(import) <- names(fit$model[2:nvar])
    names(import) <- "Weights"
    import <- import[order(import), 1, drop = FALSE]
    dotchart(import$Weights, labels = row.names(import), xlab = "% of R-Square", 
        pch = 19, main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=", 
            round(rsquare, digits = 3)), ...)
    return(import)
}
## ============================================================================##

13.2.2.3 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών

Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:

## ============================================================================##
## Check the relative importance
## ============================================================================##
rel.1 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)

relweights(rel.1, col = "blue")

	Weights
MaxD	4.671829
bio8	4.868050
monthCountByTemp10	5.001657
MinD	22.456880
bio14	31.196091
Area	31.805492

13.2.3 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων

Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:

## ============================================================================##
## The best subsets method
## ============================================================================##
bes.tot <- leaps::regsubsets(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + 
    MinD, data = ryukyus_log, nbest = 4)
## ============================================================================##

13.2.3.1 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος

Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:

## ============================================================================##
## Visualise
## ============================================================================##
plot(bes.tot, scale = "adjr2")

## ============================================================================##
## Check the summary and the relative importance
## ============================================================================##
tot <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD, data = ryukyus_log)

summary(tot)

## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89706 -0.24475 -0.02862  0.37684  0.80405 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -21.6286    11.5101  -1.879   0.0897 .
## Area                 0.7130     0.2705   2.636   0.0249 *
## bio14                7.4248     3.1244   2.376   0.0389 *
## bio8                 5.3227     3.9457   1.349   0.2071  
## monthCountByTemp10  -5.6306     5.0687  -1.111   0.2926  
## MinD                -2.1181     0.8620  -2.457   0.0338 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared:  0.6475, Adjusted R-squared:  0.4712 
## F-statistic: 3.673 on 5 and 10 DF,  p-value: 0.03798

relweights(tot, col = "blue")

	Weights
bio8	5.224379
monthCountByTemp10	5.706259
MinD	24.302432
Area	31.886016
bio14	32.880914

13.2.4 Έλεγχος των αποτελεσμάτων της παλινδρόμησης

Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:

sqrt(car::vif(tot)) > 2

##               Area              bio14               bio8 monthCountByTemp10 
##              FALSE              FALSE              FALSE              FALSE 
##               MinD 
##              FALSE

13.2.5 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

## ============================================================================##
## Compare the two models
## ============================================================================##
AIC(enter2, stepm2, tot)

	df	AIC
enter2	8	35.00027
stepm2	5	32.37180
tot	7	33.04070

13.2.6 Έλεγχος πιθανοτήτων

Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:

## ============================================================================##
## Let's have a look on the p-values for each of the IVs that have entered our
## final model
## ============================================================================##
summary(tot)

## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89706 -0.24475 -0.02862  0.37684  0.80405 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -21.6286    11.5101  -1.879   0.0897 .
## Area                 0.7130     0.2705   2.636   0.0249 *
## bio14                7.4248     3.1244   2.376   0.0389 *
## bio8                 5.3227     3.9457   1.349   0.2071  
## monthCountByTemp10  -5.6306     5.0687  -1.111   0.2926  
## MinD                -2.1181     0.8620  -2.457   0.0338 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared:  0.6475, Adjusted R-squared:  0.4712 
## F-statistic: 3.673 on 5 and 10 DF,  p-value: 0.03798

## ============================================================================##
## Inside the parenthesis we include the above p-values
## ============================================================================##
p <- c(0.0249, 0.0389, 0.2071, 0.2926, 0.0338)
## ============================================================================##


## ============================================================================##
## Level of significance
## ============================================================================##
alp <- 0.05
## ============================================================================##


## ============================================================================##
## If the model's p-values are equal or smaller than the p-values we obtain from
## the following command, then they are ok
## ============================================================================##
p.adjust(p, method = "bonferroni", n = length(p))

## [1] 0.1245 0.1945 1.0000 1.0000 0.1690

14 Θερμά σημεία βιοποικιλότητας

## ============================================================================##
## Load the biotic cleaned data, as well as the spatial data for Italy
## ============================================================================##
data_it_cleaned <- readRDS("Italian data cleaned.rds")

Italy <- read_sf("Italin provinces.shp") %>% st_transform(4326) %>% as_Spatial()

## ============================================================================##
## First load the excel file with the correct species names
## ============================================================================##
nms <- readxl::read_excel("italian_names.xlsx")
## ============================================================================##


## ============================================================================##
## Keep only one instance per taxon from our occurrence data
## ============================================================================##
nms_flags <- data_it_cleaned %>% distinct(scientificName) %>% mutate(Taxon = nms$scientificName, 
    action = nms$action) %>% as_tibble()
## ============================================================================##

## ============================================================================##
## First create the coords for the correct taxa
## ============================================================================##
coords <- data_it_cleaned %>% full_join(nms_flags) %>% as_tibble() %>% filter(!str_detect(action, 
    "DEL")) %>% dplyr::select(decimalLongitude, decimalLatitude, Taxon)
## ============================================================================##

## ============================================================================##
## Then create the sparse matrix
## ============================================================================##
pt <- points2comm(dat = coords, mask = Greece, res = 0.04166667, lon = "decimalLongitude", 
    lat = "decimalLatitude", species = "Taxon")

head(pt[[1]][1:5, 1:5])

## 5 x 5 sparse Matrix of class "dgCMatrix"
##        Abies_alba Abies_nebrodensis Abies_nordmanniana Abutilon_theophrasti
## v10000          .                 .                  .                    .
## v10003          .                 .                  .                    .
## v10004          .                 .                  .                    .
## v10007          .                 .                  .                    .
## v10008          .                 .                  .                    .
##        Acacia_dealbata
## v10000               .
## v10003               .
## v10004               .
## v10007               .
## v10008               .

## ============================================================================##
## Save it
## ============================================================================##
saveRDS(pt, "Sparse matrix and relevant polygons for Italy.rds")

##============================================================================##
## Create the phylogenetic tree
##============================================================================##
italy_tree <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(family, genus, Taxon) %>% 
  dplyr::rename(species = Taxon) %>% 
  distinct(species, .keep_all = T) %>%  
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)


##============================================================================##
## Save it to disk
##============================================================================##
saveRDS(italy_tree, 'Italian phylogenetic tree.rds')
##============================================================================##

##============================================================================##
## Locate the hotspots
##============================================================================##
hot_cold_spots <- pt$poly_shp %>% 
  
  ## Convert to sf 
  st_as_sf() %>%
  
  ## Create the hot- and cold-spots
  mutate(Cold_SR = coldspots(richness, 
                             prob = 5),
         Hot_SR = hotspots(richness, 
                           prob = 5),
         
         ## Corrected weighted endemism
         Cold_CWE = coldspots(weighted_endemism(pt$comm_dat),
                              prob = 5),
         Hot_CWE = hotspots(weighted_endemism(pt$comm_dat), 
                            prob = 5),
         
         ## Phylogenetic diversity
         Cold_PD = coldspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                             prob = 5),
         Hot_PD = hotspots(PD(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                              match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                           prob = 5),
         
         ## Phylogenetic endemism
         Cold_PE = coldspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                            match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                             prob = 5),
         Hot_PE = hotspots(phylo_endemism(match_phylo_comm(italy_tree, pt$comm_dat)$com, 
                                          match_phylo_comm(italy_tree, pt$comm_dat)$phy), 
                           prob = 5)) %>% 
  mutate(all_hot = ifelse(Hot_SR == 1 & Hot_CWE == 1 & Hot_PD == 1 & Hot_PE == 1, 'HOT', 
                          ifelse(Cold_SR == 1 & Cold_CWE == 1 & Cold_PD == 1 & Cold_PE == 1, 'COLD', 'Neither'))) %>% 
  mutate(cwe_sr_hot = ifelse(Hot_SR == 1 & Hot_CWE == 1, 'HOT', 
                          ifelse(Cold_SR == 1 & Cold_CWE == 1, 'COLD', 'Neither')))
##============================================================================##

## ============================================================================##
## Create a new variable
## ============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% as_Spatial()
## ============================================================================##


## ============================================================================##
## Quick plot
## ============================================================================##
plot(sp_hot_cold_spots, border = "grey", col = "lightgrey", main = "Diversity Hotspots and Coldspots")

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == "HOT"), ], col = "red", 
    add = TRUE, border = NA)

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == "COLD"), ], col = "blue", 
    add = TRUE, border = NA)

legend("bottomleft", fill = c("red", "blue"), legend = c("hotspots", "coldspots"), 
    bty = "n", inset = 0.092)

## ============================================================================##
## Second quick plot
## ============================================================================##
hot_cold_spots %>% ggplot() + geom_sf(aes(fill = cwe_sr_hot)) + labs(title = "Diversity hotspots in Italy", 
    subtitle = "Source: GBIF") + scale_fill_manual(values = c("steelblue1", "red", 
    "gray90"), name = "Diversity hotspots", labels = c("COLD", "HOT", "Neither"))

15 Appendix: R code

##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
options(Encoding="UTF-8")

library(knitr)
library(rmdformats)

## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
                 cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               # comment=NA,
               message=FALSE,
               warning=FALSE,
               eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
htmltools::img(src = knitr::image_uri(file.path("G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/HBS Seminar/Logo_HBS/ebe-LOGO-005-300_1_70.jpg")), 
               alt = 'logo', 
               style = 'position:absolute; top:0; right:0; padding:5px;')
##===========================================================================##
## Load the libraries
##===========================================================================##
library(tidyverse)
library(raster)
library(sf)
library(ape)
library(phyloregion)
library(picante)
library(rtrees)
library(ggtree)
library(exactextractr)
library(caret)
library(psych)
library(PhyloMeasures)
library(phyloregion)
library(Matrix)
library(broom)
##===========================================================================##
##===========================================================================##
## Load our data
##===========================================================================##
plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", row.names = 1)
##===========================================================================##
##===========================================================================##
## Check their dimensions
##===========================================================================##
dim(plants) 
##===========================================================================##
## Check their dimensions
##===========================================================================##
str(plants) 
pre {
  max-height: 400px;
  overflow-y: auto;
}

pre[class] {
  max-height: 400px;
}
##===========================================================================##
## Convert the P/A matrix to SR
##===========================================================================##
ryukyus <- rowSums(plants) %>% 
  enframe() %>% 
  dplyr::rename(Island = name, 
                SR = value)

ryukyus
##===========================================================================##
## Create a vector containing just the file paths
##===========================================================================##
shps <- fs::dir_ls("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/", 
                   glob = "*.shp")
##===========================================================================##


##===========================================================================##
## Create a vector containing just the island names
##===========================================================================##
island_nms <- shps %>% 
  str_replace_all('C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/',
                  '') %>% 
  str_replace_all('\\.shp', '')
##===========================================================================##


##===========================================================================##
## Bulk load the shapefiles
##===========================================================================##
ryukyus_shps <- shps %>% 
  purrr::map(read_sf) %>% 
  lapply(., dplyr::select, NAME_1, geometry)  %>% 
  do.call(rbind, .) %>% 
  mutate(Island = island_nms) %>% 
  dplyr::select(Island, geometry)
##===========================================================================##
##===========================================================================##
## Load, crop and mask the rasters
##===========================================================================##
preds <- stack(stack(list.files(path = "I:/SDM DATA/WORLDCLIM/Current/bio_2-5m_bil/",
                                pattern = '.bil',
                                recursive = T,
                                full.names = T)) %>% 
                 crop(., ryukyus_shps, snap = 'in') %>% 
                 mask(., ryukyus_shps),
               stack(list.files(path = "I:/SDM DATA/Envirem_Predictors/Eurasia_current_2.5arcmin_geotiff/",
                                pattern = '.tif',
                                recursive = T,
                                full.names = T))  %>% 
                 crop(., ryukyus_shps, snap = 'in') %>% 
                 mask(., ryukyus_shps))
##===========================================================================##
##===========================================================================##
## Load, crop and mask the altitudinal raster
##===========================================================================##
altitude <- raster("I:/SDM DATA/Altitude/elevation_1KMmd_GMTEDmd.tif")  %>% 
  crop(., ryukyus_shps, snap = 'in') %>% 
  mask(., ryukyus_shps)
##===========================================================================##
##===========================================================================##
## Extract abiotic info
##===========================================================================##
preds_ryukyus <- exact_extract(preds, ryukyus_shps, 'median') %>% 
  mutate(Island = ryukyus_shps$Island) %>% 
  full_join(.,
            exact_extract(altitude, ryukyus_shps, 'median') %>% 
              enframe() %>% 
              mutate(Island = ryukyus_shps$Island) %>% 
              rename(Elevation = value) %>% 
              dplyr::select(Island, Elevation)) %>% 
  dplyr::select(Island, everything()) %>% 
  as_tibble()


##============================================================================##
## Give them proper names
##============================================================================##
col_names <- preds_ryukyus %>%
  colnames() %>%
  enframe() %>%
  mutate(value = str_replace_all(value, 'median\\.', '')) %>%
  mutate(value = str_replace_all(value, 'current_2\\.5arcmin_', '')) %>% 
  pull(value)
##===========================================================================##


##============================================================================##
## Overwrite names and remove Row 1
##============================================================================##
vars <- preds_ryukyus %>%
  set_names(col_names) %>% 
  dplyr::select(Island, everything()) %>% 
  as_tibble()
##============================================================================##
## Transform the missing data
##============================================================================##
imputed_vars <- predict(preProcess(as.data.frame(vars[, 2:36]), 
                                   method = c('medianImpute')), 
                        vars)
##============================================================================##
## First drop any rows with NA values
##============================================================================##
drop_vars <- vars %>% drop_na()
##============================================================================##


##============================================================================##
## Then check which islands have NA values in any of the variables
##============================================================================##
setdiff(imputed_vars$Island, drop_vars$Island)


##============================================================================##
## Create a symbol to filter out all the Islands having NA values
##============================================================================##
'%!in%' <- function(x,y)!('%in%'(x,y))
##============================================================================##


##============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
##============================================================================##
plants_new <- plants %>% 
  mutate(Island = rownames(.)) %>% 
  as_tibble() %>% 
  dplyr::select(Island, everything()) %>% 
  pivot_longer(!Island, names_to = 'Taxon', 
               values_to = 'PA') %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  pivot_wider(names_from = Taxon,
              values_from = PA)
##============================================================================##


##============================================================================##
## See which taxa are left after removing the problematic islands
##============================================================================##
keepers <- colSums(plants_new[, -1]) %>% 
  enframe() %>% 
  filter(value > 0)
##============================================================================##


##============================================================================##
## Finally create the new PA matrix
##============================================================================##
plants_new_pa <- plants_new %>% 
  pivot_longer(!Island, names_to = 'Taxon', 
               values_to = 'PA') %>% 
  filter(Taxon %in% keepers$name) %>% 
  pivot_wider(names_from = Taxon,
              values_from = PA)
##============================================================================##
datatable(plants, options = list(pagelength = 10)) 
jpn_tree <- readRDS('G:/Academic/Lessons/EMB/Labs/Labs/16 Ryukyus tree.rds')
##============================================================================##
## Plot the tree
##============================================================================##
groupInfo <- split(jpn_tree$tip.label, 
                   gsub("_\\w+", "", 
                        jpn_tree$tip.label))

groups_jpn <- groupOTU(jpn_tree, 
                       groupInfo)
ggtree(groups_jpn, aes(color = group),
       layout = 'circular') +
  geom_tiplab2(size = 0.5, align = T) 
##============================================================================##
## Combine the trees with the P/A data
##============================================================================##

#check for mismatches/missing species
combined_japan <- picante::match.phylo.comm(jpn_tree,plants_new_pa[, 2:1800])

##============================================================================##
## The resulting object is a list with $phy and $comm elements
## replace our original data with the sorted/matched data
##============================================================================##
phy_japan <- combined_japan$phy

comm_japan <- combined_japan$comm
##============================================================================##
##============================================================================##
## Compute phylogenetic alpha diversity
##============================================================================##

## For the estimation of SES values
japan_pd <- pd.query(phy_japan, comm_japan, TRUE) 

## For the p-values
japan_pd_p <- pd.pvalues(phy_japan, comm_japan, reps = 1000)
##============================================================================##
##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>% 
  filter(Island %in% plants_new_pa$Island) %>% 
  full_join(., drop_vars) %>% 
  full_join(., ryukyus_shps) %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  st_as_sf() %>% 
  
  ## Lambert Azimuthal Equal Area projection
  st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>% 
  
  ## Create the Area and PD variables
  mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
         PD = japan_pd,
         PD_P = japan_pd_p) %>% 
  
  ## Transform back to WGS84
  st_transform(4326) %>% 
  
  dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>% 
  
  ## It has zero-variance
  dplyr::select(-aridityIndexThornthwaite) %>% 
  
  ## Drop the geometry
  st_drop_geometry() %>% 
  
  ## Join these data with the Maximum and Minimum distance between each island
  ## pair
  inner_join(., 
             
             ## Calculate the Maximum distance between each island pair
             apply((st_distance(ryukyus_shps %>% 
                                     st_centroid(), 
                                   ryukyus_shps %>% 
                                     st_centroid())/1e+3) %>% 
                        as_tibble() %>% 
                        set_names(ryukyus_shps$Island) %>% 
                        units::drop_units() %>% 
                        na_if(., 0), 1, max, na.rm = T) %>% 
               enframe() %>% 
               rename(MaxD = value) %>% 
               mutate(Island = ryukyus_shps$Island) %>% 
               dplyr::select(Island, MaxD) %>% 
               
               full_join(., 
                         
                         ## Calculate the Minimum distance between each island
                         ## pair
                         apply((st_distance(ryukyus_shps %>% 
                                                 st_centroid(), 
                                               ryukyus_shps %>% 
                                                 st_centroid())/1e+3) %>% 
                                    as_tibble() %>% 
                                    set_names(ryukyus_shps$Island) %>% 
                                    units::drop_units() %>% 
                                    na_if(., 0), 1, min, na.rm = T) %>% 
                           enframe() %>% 
                           rename(MinD = value) %>% 
                           mutate(Island = ryukyus_shps$Island) %>% 
                           dplyr::select(Island, MinD)))
##============================================================================##

ryukyus_new
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- keepers %>% 
  separate(name, into = c('genus', 'binomial'), sep = '_', remove = F) %>% 
  rename(species = name) %>% 
  dplyr::select(species, genus) %>% 
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)
##============================================================================##
##============================================================================##
## Normality check
##============================================================================##
normality <- lapply(ryukyus_new[,2:42], shapiro.test)
##============================================================================##
normality
##============================================================================##
## Correlation check
##============================================================================##
corr.test(ryukyus_new[,5:42], method = "spearman")
##============================================================================##
## Visualization
##============================================================================##
car::scatterplotMatrix(ryukyus_new[,5:10], 
                  spread = F,
                  smoother.args = list(lty = 2), 
                  main = "Scatter Plot Matrix")
##============================================================================##
## Correlation check
##============================================================================##
unconc <- usdm::vifcor(ryukyus_new %>% 
                         as.data.frame() %>% 
                         dplyr::select(Area:MinD), 
                       th = 0.7)
unconc <- unconc@results$Variables[unconc@results$VIF <= 5]
unconc <- droplevels(factor(unconc))
##============================================================================##
##============================================================================##
## Log-transform the data
##============================================================================##
ryukyus_log <- ryukyus_new %>% 
  dplyr::select(Area, bio14, bio8,
                monthCountByTemp10, 
                MaxD, MinD) %>% 
  mutate(SR = ryukyus_new$SR) %>% 
  log10(.) %>% 
  mutate(Island = ryukyus_new$Island,
         PD = ryukyus_new$PD) %>% 
  dplyr::select(Island, SR, PD, everything())
##============================================================================##
##============================================================================##
## ISAR models
##============================================================================##
model1 <- lm(SR ~ Area, data = ryukyus_log)
model2 <- lm(PD ~ Area, data = ryukyus_log)
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(model1)
##============================================================================##
## Check the summary
##============================================================================##
summary(model2)
##============================================================================##
## Create the residuals
##============================================================================##
res_sr <- rstandard(model1)
res_pd <- rstandard(model2)
##============================================================================##
##============================================================================##
## Check the estimate and the p-value
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% tidy()) %>% 
  mutate(Metric = rep(c('SR','PD'), each = 2)) %>% 
  dplyr::select(Metric, everything())
##============================================================================##

##============================================================================##
## Check the metrics
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% glance()) %>% 
  mutate(Metric = c('SR','PD')) %>% 
  dplyr::select(Metric, everything())
##============================================================================##


##============================================================================##
## Check the residuals
##============================================================================##
resids <- lm(SR ~ Area, 
             data = ryukyus_log) %>% 
  augment() %>% 
  mutate(Island = ryukyus_log$Island)  %>% 
  rename(SR_res = .std.resid) %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% 
              augment() %>% 
              mutate(Island = ryukyus_log$Island) %>% 
              rename(PD_res = .std.resid),
            by = 'Island') %>% 
  dplyr::select(Island, SR, SR_res, PD, PD_res)
##============================================================================##
##============================================================================##
## Plot
##============================================================================##
ggplot(ryukyus_log, aes(res_sr, res_pd)) +
  geom_point() + 
  geom_text(aes(label = Island), 
            fontface = "bold", 
            size = 2.8, 
            vjust = -1.2, 
            hjust = 0.5) + 
  geom_hline(yintercept = 0.0) + 
  geom_vline(xintercept = 0.0) + 
  theme(legend.position = "none") + 
  xlab("SR") + 
  ylab("PD") 
##============================================================================##
## Quick save
##============================================================================##
ggsave("SR-PD.png", width = 20, height = 20, units = "cm", dpi = 600)
##============================================================================##
## The enter method
##============================================================================##
enter1 <- lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log) 

enter2 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log) 
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(enter1)
summary(enter2)

##============================================================================##
## The solution
##============================================================================##
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm2)
##============================================================================##
## The stepAIC method
##============================================================================##
stepm1 <- MASS::stepAIC(lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log), 
                  scale = 0, 
                  direction = "both", 
                  trace = 1, 
                  keep = NULL, 
                  steps = 1000,
                  k = 2)

stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm1)
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm2)
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter1, stepm1)
##============================================================================##
## A function for the relative importance
##============================================================================##
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import),1, drop=FALSE]
dotchart(import$Weights, labels=row.names(import),
xlab="% of R-Square", pch=19,
main="Relative Importance of Predictor Variables",
sub=paste("Total R-Square=", round(rsquare, digits=3)),
...)
return(import)
}
##============================================================================##
##============================================================================##
## Check the relative importance
##============================================================================##
rel.1 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)

relweights(rel.1, col="blue")
##============================================================================##
## The best subsets method
##============================================================================##
bes.tot <- leaps::regsubsets(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, 
                      data = ryukyus_log, 
                      nbest = 4)
##============================================================================##
##============================================================================##
## Visualise
##============================================================================##
plot(bes.tot, scale = "adjr2")
##============================================================================##
## Check the summary and the relative importance
##============================================================================##
tot <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD, data = ryukyus_log)  

summary(tot)


relweights(tot, col="blue")
sqrt(car::vif(tot)) > 2
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter2, stepm2, tot)
##============================================================================##
## Load the biotic cleaned data, as well as the spatial data for Italy                            
##============================================================================##
data_it_cleaned <- readRDS('Italian data cleaned.rds')

Italy <- read_sf('Italin provinces.shp') %>% 
  st_transform(4326) %>% 
  as_Spatial()
##============================================================================##
## First load the excel file with the correct species names
##============================================================================##
nms <- readxl::read_excel('italian_names.xlsx')
##============================================================================##


##============================================================================##
## Keep only one instance per taxon from our occurrence data
##============================================================================##
nms_flags <- data_it_cleaned %>% 
  distinct(scientificName) %>% 
  mutate(Taxon = nms$scientificName,
         action = nms$action) %>% 
  as_tibble() 
##============================================================================##

##============================================================================##
## First create the coords for the correct taxa
##============================================================================##
coords <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(decimalLongitude, 
                decimalLatitude,
                Taxon)
##============================================================================##


pt <- readRDS('Sparse matrix and relevant polygons for Italy.rds')
##============================================================================##
## Then create the sparse matrix
##============================================================================##
pt <- points2comm(dat = coords, 
                  mask = Greece, 
                  res = 0.04166667, 
                  lon = "decimalLongitude", 
                  lat = "decimalLatitude",
                  species = "Taxon")
head(pt[[1]][1:5, 1:5])
##============================================================================##
## Save it
##============================================================================##
saveRDS(pt, 'Sparse matrix and relevant polygons for Italy.rds')
italy_tree <- readRDS('Italian phylogenetic tree.rds')
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
italy_tree <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(family, genus, Taxon) %>% 
  dplyr::rename(species = Taxon) %>% 
  distinct(species, .keep_all = T) %>%  
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)


##============================================================================##
## Save it to disk
##============================================================================##
saveRDS(italy_tree, 'Italian phylogenetic tree.rds')
##============================================================================##
##============================================================================##
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##============================================================================##
summary(tot) 



##============================================================================##
## Inside the parenthesis we include the above p-values             
##============================================================================##
p <- c(0.0249, 0.0389, 0.2071, 0.2926, 0.0338) 
##============================================================================##


##============================================================================##
## Level of significance                             
##============================================================================##
alp <- 0.05 
##============================================================================##


##============================================================================##
## If the model's p-values are equal or smaller than the p-values we obtain 
## from the following command, then they are ok
##============================================================================##
p.adjust(p, method = "bonferroni", n = length(p)) 
##============================================================================##
## Create a new variable
##============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% 
  as_Spatial()
##============================================================================##


##============================================================================##
## Quick plot
##============================================================================##
plot(sp_hot_cold_spots,
     border = "grey",
     col = "lightgrey",
     main = "Diversity Hotspots and Coldspots")

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'HOT'), ],
     col = "red",
     add = TRUE, 
     border = NA)

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'COLD'), ], 
     col = "blue", 
     add = TRUE, 
     border = NA)

legend("bottomleft",
       fill = c("red", "blue"),
       legend = c("hotspots", "coldspots"),
       bty = "n", 
       inset = .092)


##============================================================================##
## Second quick plot
##============================================================================##
hot_cold_spots %>% 
  ggplot() +
  geom_sf(aes(fill = cwe_sr_hot )) + 
  labs(title = 'Diversity hotspots in Italy',
       subtitle = 'Source: GBIF') + 
  scale_fill_manual(values = c('steelblue1','red', 'gray90'),
                    name = "Diversity hotspots",
                    labels = c("COLD", "HOT", "Neither"))
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##